!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!!
!! Semidiscretized ODE for the compressible Navier-Stokes DG system
!!
!! \n
!!
!! This module represents the semi-discretized problem as an ODE, with
!! the layout required by the time integrators of \c
!! mod_time_integrators.
!!
!! \section turbulence_model Turbulence model
!!
!! In general, a complex turbulence model requires:
!! <ul>
!!  <li> additional terms in the right-hand-side of the NS equations
!!  <li> additional diagnostic quantities, for instance the gradients
!!  of some prognostic variables
!!  <li> additional prognostic variables, such as time averages of
!!  other quantities or variables related to the turbulent kinetic
!!  energy for \f$\kappa-\epsilon\f$ and \f$\kappa-\omega\f$ models.
!! </ul>
!! These objects have to be included in the general layout for ODE
!! problems provided by \c mod_time_integrators. We proceed as
!! follows.
!! <ul>
!!  <li> The "model" is stored in <tt>t_dgcomp_ode\%turbmod</tt>:
!!  this includes a collection of subroutines to compute the
!!  additional terms in the equations, as well as the tendencies of
!!  the additional prognostic variables in case they are present;
!!  these subroutines can be called in \c mod_dgcomp_rhs (notice that
!!  the \c t_dgcomp_ode object is passed to the computation of the
!!  rhs). <tt>t_dgcomp_ode\%turbmod</tt> also includes subroutines for
!!  initialization and clean-up of the turbulence model related
!!  objects included in \c c_stv and \c c_ods.
!!  <li> The additional diagnostic quantities are included in
!!  <tt>t_dgcomp_ods\%turbmod_diags</tt>. Since these diagnostics are
!!  model specific, they must be defined by a dedicated method of
!!  <tt>t_dgcomp_ode\%turbmod</tt>. Notice that the creation of new \c
!!  \c t_dgcomp_stv, \c t_dgcomp_ode and \c t_dgcomp_ods objects is
!!  not considered in this module, except for some minimal
!!  initialization, so all these initializations should be done
!!  elsewhere.
!!  <li> The additional prognostic variables are included in
!!  <tt>t_dgcomp_stv\%turbmod_progs</tt>. The treatment of this object
!!  is very similar to <tt>t_dgcomp_ods\%turbmod_diags</tt>.
!! </ul>
!<----------------------------------------------------------------------
module mod_dgcomp_ode

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized, &
   wp_mpi, mpi_sum, mpi_allgather

 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

 use mod_time_integrators, only: &
   mod_time_integrators_initialized, &
   c_ode, c_ods

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid

 use mod_bcs, only: &
   mod_bcs_initialized, &
   t_bcs

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_turb_flux, only: &
   mod_turb_flux_initialized, &
   c_turbmod, c_turbmod_progs, c_turbmod_diags, t_turb_diags

 use mod_dgcomp_rhs, only: &
   mod_dgcomp_rhs_initialized, &
   t_bcs_error, &
   dgcomp_tens

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_dgcomp_ode_constructor, &
   mod_dgcomp_ode_destructor,  &
   mod_dgcomp_ode_initialized, &
   t_dgcomp_stv, t_dgcomp_ods, t_dgcomp_ode, &
   new_dgcomp_stv, new_dgcomp_ode, clear

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> State information for the compressible Navier-Stokes equations.
 !!
 !! Some pointers to the finite element space description are also
 !! introduced; such description is shared by all the state
 !! variables.
 type, extends(c_stv) :: t_dgcomp_stv
  type(t_grid), pointer :: grid => null()
  type(t_base), pointer :: base => null()
  type(t_bcs ), pointer :: bcs  => null()
  !> The state vector is organized on three indexes: the first one
  !! identifies the variable, the second one the element degree of
  !! freedom, the third on the mesh element
  real(wp), allocatable :: u(:,:,:)
  !> Integrated mass flow
  !! \f{displaymath}{
  !!  \underline{\mathcal{M}}_{\mathcal{I}} = \int_{t_0}^t \left(
  !!  \underline{M}(s) - \bar{\underline{M}}\right)ds
  !! \f}
  !! where the mass flow \f$\underline{M}\f$ is defined in \c
  !! mod_dgcomp_flowstate and \f$\bar{\underline{M}}\f$ is a reference
  !! value specified by the selected test case. Clearly, the
  !! corresponding evolution equation is
  !! \f{displaymath}{
  !!  \dot{\underline{\mathcal{M}}}_{\mathcal{I}} = 
  !!  \underline{M} - \bar{\underline{M}}.
  !! \f}
  !! \note The mass flow \f$\underline{M}\f$ is integrated over the
  !! whole volume, and so is
  !! \f$\underline{\mathcal{M}}_{\mathcal{I}}\f$. This means that all
  !! processors share a copy of the same value.
  real(wp), allocatable :: im(:)
  !> This field can be used to include additional equations in the
  !! problem, for instance turbulence models.
  !!
  !! This field can be left unallocated, in which case it is ignored
  !! in this module. If it is allocated, then the corresponding
  !! operations are used: for instance, <tt>t_dgcomp_stv\%incr</tt>
  !! calls <tt>t_dgcomp_stv\%turbmod_progs\%incr</tt> and so on.
  !!
  !! This field is not allocated here: if you want to use it, it must
  !! be allocated after defining a new \c t_dgcomp_stv object with \c
  !! new_dgcomp_stv.
  class(c_turbmod_progs), allocatable :: turbmod_progs
  integer :: comm
  type(t_bcs_error) :: bcerr
 contains
  procedure, pass(x) :: incr
  procedure, pass(x) :: tims
  !> This subroutine is overidden forbetter efficiency
  procedure, pass(x) :: inlt => dgcomp_inlt
  procedure, pass(z) :: copy
  procedure, pass(x) :: scal => dgcomp_scal
  procedure, pass(x) :: source
  procedure, pass(x) :: source_vect
 end type t_dgcomp_stv

 !> ODE for the compressible Navier-Stokes equations.
 type, extends(c_ode) :: t_dgcomp_ode
  type(t_grid), pointer :: grid => null()
  type(t_base), pointer :: base => null()
  type(t_bcs ), pointer :: bcs  => null()
  logical :: viscous_flow
  integer :: comm
  class(c_turbmod), allocatable :: turbmod
 contains
  procedure, pass(ode) :: rhs  => dgcomp_ode_rhs
 end type t_dgcomp_ode

 !> ODE diagnostic variables
 type, extends(c_ods) :: t_dgcomp_ods
  class(c_turbmod_diags), allocatable :: turbmod_diags
 end type t_dgcomp_ods
 
 ! private members

 interface clear
   module procedure clear_dgcomp_stv, clear_dgcomp_ode
 end interface

! Module variables
 real(wp) :: sendbuf(1)
 real(wp), allocatable :: recvbuf(:)

 ! public members
 logical, protected ::               &
   mod_dgcomp_ode_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_dgcomp_ode'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_dgcomp_ode_constructor(mpi_nd)
  integer, intent(in) :: mpi_nd

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
   (mod_mpi_utils_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
  (mod_state_vars_initialized.eqv..false.) .or. &
(mod_time_integrators_initialized.eqv..false.) .or. &
        (mod_grid_initialized.eqv..false.) .or. &
         (mod_bcs_initialized.eqv..false.) .or. &
        (mod_base_initialized.eqv..false.) .or. &
   (mod_turb_flux_initialized.eqv..false.) .or. &
  (mod_dgcomp_rhs_initialized.eqv..false.) ) then 
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_dgcomp_ode_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   allocate(recvbuf(mpi_nd))

   mod_dgcomp_ode_initialized = .true.
 end subroutine mod_dgcomp_ode_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_dgcomp_ode_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_dgcomp_ode_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   deallocate(recvbuf)

   mod_dgcomp_ode_initialized = .false.
 end subroutine mod_dgcomp_ode_destructor

!-----------------------------------------------------------------------
 
 subroutine new_dgcomp_stv(obj,grid,base,bcs,nvar,comm)
  type(t_grid), intent(in), target :: grid
  type(t_base), intent(in), target :: base
  type(t_bcs ), intent(in), target :: bcs
  integer, intent(in) :: nvar
  integer, intent(in) :: comm
  type(t_dgcomp_stv), intent(out) :: obj

   obj%grid => grid
   obj%base => base
   obj%bcs  => bcs
   allocate( obj%u(nvar,base%pk,grid%ne) )
   allocate( obj%im(grid%d) )

   ! the field obj%turbmod_progs is not allocated here; if you need it,
   ! you should allocate it after calling this subroutine

   obj%comm = comm

   ! leaving the default initialization of bcerr

 end subroutine new_dgcomp_stv
 
!-----------------------------------------------------------------------
 
 pure subroutine clear_dgcomp_stv(obj)
  type(t_dgcomp_stv), intent(out) :: obj

   nullify( obj%grid )
   nullify( obj%base )
   nullify( obj%bcs  )
   ! allocatable components are implicitly deallocated
 
 end subroutine clear_dgcomp_stv
 
!-----------------------------------------------------------------------

 subroutine incr(x,y)
  class(c_stv),        intent(in)    :: y
  class(t_dgcomp_stv), intent(inout) :: x

   select type(y); type is(t_dgcomp_stv)

   x%u = x%u + y%u
   x%im= x%im+ y%im

   if(allocated(x%turbmod_progs)) &
     call x%turbmod_progs%incr(y%turbmod_progs)

   ! propagate bcs errors (errors on x have larger priority)
   if( (.not.x%bcerr%lerr) .and. y%bcerr%lerr ) &
     x%bcerr = y%bcerr

   end select
 end subroutine incr

!-----------------------------------------------------------------------

 subroutine tims(x,r)
  real(wp),            intent(in)    :: r
  class(t_dgcomp_stv), intent(inout) :: x

   x%u = r*x%u
   x%im= r*x%im
   if(allocated(x%turbmod_progs)) call x%turbmod_progs%tims(r)

 end subroutine tims

!-----------------------------------------------------------------------

 subroutine dgcomp_inlt(x,r,y)
  real(wp),     intent(in) :: r
  class(c_stv), intent(in) :: y
  class(t_dgcomp_stv), intent(inout) :: x

   select type(y); type is(t_dgcomp_stv)

   x%u = x%u + r*y%u
   x%im= x%im+ r*y%im
   if(allocated(x%turbmod_progs)) &
     call x%turbmod_progs%inlt(r,y%turbmod_progs)

   ! propagate bcs errors (errors on x have larger priority)
   if( (.not.x%bcerr%lerr) .and. y%bcerr%lerr ) &
     x%bcerr = y%bcerr

   end select

 end subroutine dgcomp_inlt

!-----------------------------------------------------------------------

 subroutine copy(z,x)
  class(c_stv),        intent(in)    :: x
  class(t_dgcomp_stv), intent(inout) :: z

   select type(x); type is(t_dgcomp_stv)
   z%u     = x%u
   z%im    = x%im
   if(allocated(z%turbmod_progs)) &
     call z%turbmod_progs%copy(x%turbmod_progs)
   z%bcerr = x%bcerr
   end select
 end subroutine copy

!-----------------------------------------------------------------------

 function dgcomp_scal(x,y) result(s)
  class(t_dgcomp_stv), intent(in) :: x
  class(c_stv),        intent(in) :: y
  real(wp) :: s

  integer :: ierr
  character(len=*), parameter :: &
    this_fun_name = 'dgcomp_ode_scal'

   select type(y); type is(t_dgcomp_stv)

   ! This is easy because all dofs are local
   sendbuf(1) = sum(x%u*y%u)
   ! The proper form should use mpi_allreduce, however this has
   ! problem with quadruple precision:
   !   https://trac.mcs.anl.gov/projects/mpich2/ticket/258
   ! Thus we communicate all the values and use the fortran intrinsic
   ! sum to compute the reduction
   !call mpi_allreduce(sendbuf,recvbuf,1,wp_mpi,mpi_sum,x%comm,ierr)
   call mpi_allgather(sendbuf,1,wp_mpi,recvbuf,1,wp_mpi,x%comm,ierr)
   s = sum(recvbuf)

   ! Integrated mass flow: as discussed in mod_dgcomp_ode, all
   ! processors share a copy of the same value, so this contribution
   ! can be added locally
   s = s + sum( x%im * y%im )

   if(allocated(x%turbmod_progs)) &
     s = s + x%turbmod_progs%scal(y%turbmod_progs)

   if( x%bcerr%lerr .or. y%bcerr%lerr ) s = -huge(s)

   end select
 end function dgcomp_scal

!-----------------------------------------------------------------------

 subroutine source(y,x)
  class(t_dgcomp_stv), intent(in) :: x
  class(c_stv), allocatable, intent(out) :: y

   allocate(t_dgcomp_stv::y)
 
   select type(y); type is(t_dgcomp_stv)

   call new_dgcomp_stv( y , x%grid,x%base,x%bcs , &
                             size(x%u,1) , x%comm )
   if(allocated(x%turbmod_progs)) &
     allocate(y%turbmod_progs,source=x%turbmod_progs)

   end select
 end subroutine source
 subroutine source_vect(y,x,m)
  integer, intent(in) :: m
  class(t_dgcomp_stv), intent(in)  :: x
  class(c_stv), allocatable, intent(out) :: y(:)

  integer :: i

   allocate(t_dgcomp_stv::y(m))

   select type(y); type is(t_dgcomp_stv)

   do i=1,m
     call new_dgcomp_stv( y(i) , x%grid,x%base,x%bcs , &
                                  size(x%u,1) , x%comm )
     if(allocated(x%turbmod_progs)) &
       allocate(y(i)%turbmod_progs,source=x%turbmod_progs)
   enddo

   end select
 end subroutine source_vect

!-----------------------------------------------------------------------

 subroutine new_dgcomp_ode(obj,grid,base,bcs,viscous_flow,comm)
  type(t_grid), intent(in), target :: grid
  type(t_base), intent(in), target :: base
  type(t_bcs ), intent(in), target :: bcs
  logical, intent(in) :: viscous_flow
  integer, intent(in) :: comm
  type(t_dgcomp_ode), intent(out) :: obj

   obj%grid => grid
   obj%base => base
   obj%bcs  => bcs
   obj%viscous_flow = viscous_flow
   obj%comm = comm

   ! the field obj%turbmod is not allocated here; if you need it, you
   ! should allocate it after calling this subroutine

 end subroutine new_dgcomp_ode
 
!-----------------------------------------------------------------------

 pure subroutine clear_dgcomp_ode(obj)
  type(t_dgcomp_ode), intent(out) :: obj

   nullify( obj%grid )
   nullify( obj%base )
   nullify( obj%bcs  )
   ! allocatable components are implicitly deallocated
 
 end subroutine clear_dgcomp_ode
 
!-----------------------------------------------------------------------


 subroutine dgcomp_ode_rhs(tnd,ode,t,uuu,ods,term)
  class(t_dgcomp_ode), intent(in) :: ode
  real(wp),     intent(in)    :: t
  class(c_stv), intent(in)    :: uuu
  class(c_ods), intent(inout) :: ods
  class(c_stv), intent(inout) :: tnd
  integer,      intent(in), optional :: term(2)

   select type(uuu); type is(t_dgcomp_stv)
     select type(tnd); type is(t_dgcomp_stv)
       select type(ods); type is(t_dgcomp_ods)

   call dgcomp_tens(tnd%u, tnd%im, tnd%turbmod_progs, tnd%bcerr, &
          ode%grid, ode%base, ode%bcs, ode%viscous_flow, ode%turbmod, &
                 t, uuu%u, uuu%im, uuu%turbmod_progs, &
                    ods%turbmod_diags, term=term)

       end select
     end select
   end select
 end subroutine dgcomp_ode_rhs

!-----------------------------------------------------------------------

end module mod_dgcomp_ode

